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Abstract 

We critically compare the practicality and accuracy of numerical approximations 
of phase field models and sharp interface models of solidification. Particular empha- 
sis is put on Stefan problems, and their quasi-static variants, with applications to 
crystal growth. New approaches with a high mesh quality for the parametric approx- 
imations of the resulting free boundary problems and new stable discretizations of 
the anisotropic phase field system are taken into account in a comparison involving 
benchmark problems based on exact solutions of the free boundary problem. 
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1 Introduction 

The solidification of a liquid or the melting of a solid lead to complex free boundary prob- 
lems involving many different physical effects. For example, latent heat is set free at the 
interface and a competition between surface energy and diffusion leads to instabilities like 
the MuUins-Sekerka instability. The resulting model is a Stefan problem with boundary 
conditions taking surface energy effects and kinetic effects at the interface into account, see 
e.g. [49l[35]. Crystals forming in an undercooled melt lead to very complex patterns and, in 
particular, dendritic growth can be observed since the growth is typically diffusion limited, 
see [19]. 
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The numerical simulation of time-dependent Stefan problems, or variants of it, is a 
formidable task since the evolving free boundary has to be computed. Hence, direct front 
tracking type numerical methods need to adequately capture the geometry of the interface 
and to evolve the interface approximation, often with a coupling to other physical fields. 
This coupling, in particular, represents a significant initial hurdle towards obtaining prac- 
tical implementations, and thus numerical simulations for the problem at hand. Examples 
of the implementation of such direct methods can be found in e.g. [SU [TTJ [5l [21 [69l [68l 1701 
EU [71102]. 

A further drawback of direct front tracking methods has been the inability of most 
direct methods to deal with so-called mesh eflFects, or to prevent them altogether. When a 
discrete approximation of an interface, for example a polygonal curve in the plane, evolves 
in time, then in general it is possible for the approximation to deteriorate or to break down. 
Examples of such pathologies are self-crossings and vertex coalescence. While for simple 
isotropic problems in the plane these issues can be dealt with, for example by frequent 
remeshings or by using clever formulations that only allow equidistributed approximations, 
see e.g. [54l IT^, until very recently there has been no remedy for fully anisotropic problems 
in two and three space dimensions. 

However, building on their work for isotropic problems in [HI [71 [10], the present authors 
recently introduced stable parametric finite element schemes for the direct approximation of 
anisotropic geometric evolution equations in [9i [11] , for which good mesh properties can be 
guaranteed. In particular, even for the simulation of interface evolutions in the presence of 
strong anisotropics, no remeshing or redistribution of vertices is needed in practice. These 
schemes, in which only the interface without a coupling to bulk quantities is modelled, 
have been extended to approximations of the Stefan problem with fully anisotropic Gibbs- 
Thomson law and kinetic undercooling in [J^. The novel method from [12] can be shown 
to be stable and to have good mesh properties. We remark that these approaches extend 
earlier ideas from [371 [69l [ZQ] . Here we recall the pioneering work of Schmidt [69l[70], where 
the full Stefan problem in three dimensions was solved within a sharp interface framework 
for the first time. 

Phase field methods are an alternative approach to solve solidification phenomena in the 
framework of continuum modelling. In phase field approaches a new non-conserved order 
parameter is introduced, which in the two phases is close to two different prescribed values 
and which smoothly changes its value across a small diffuse interfacial region. A parabolic 
partial differential equation for Lp is then coupled to an energy balance, which results in a 
diffusion equation for the temperature taking latent heat effects into account. We refer to 
[Ml [SI [251 [Ml [79] and to the five review articles [SI [SH [Ml [ZSl E5| for further details. In 
particular, it can be shown that solutions to the phase field equations converge to classical 
sharp interface problems, see e.g. [261 [H [TH [27] and the references therein. 

The popularity of phase field methods, often also called diffuse interface methods, can 
be explained by two characteristic features that they share with the level set method, which 
is another sharp interface computational tool. Firstly, phase field methods naturally allow 
for changes of topology. And secondly, computing simulations for phase field models only 
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requires the solution of partial differential equations on standard Cartesian domains. The 
fact that these can usually be implemented and solved in a relatively straightforward way 
makes the phase field method particularly appealing. 

It is the aim of this article to investigate the practical merits of phase field models 
compared to the recently introduced sharp interface algorithm for the approximation of 
Stefan problems from [12j; see also [131 El- Of particular interest will be the relative 
accuracy of the two methods, in situations where a true solution to the sharp interface 
problem is known. In a phase field simulation the observed error is made up of contributions 
from 

• the asymptotic error, 

• the spatial and temporal discretization errors, 

• rounding errors and solver tolerances. 

Here the asymptotic error is induced by the choice of interfacial parameter e > 0. In general 
one can formally show that the asymptotic error converges to zero as 6: ^ 0, see e.g. [25j. 
For certain phase field models and under certain conditions it can be rigorously shown that 
the asymptotic error vanishes as 6: ^ 0, see e.g. In computations for sharp interface 

approximations, on the other hand, the observed error is made up of the latter two contri- 
butions only, i.e. of discretization and rounding errors. A disadvantage of phase field models 
is that the resulting PDEs become stiflF for decreasing ^, leading to a requirement for very 
fine spatial and temporal discretizations. Hence it becomes computationally challenging 
to produce very accurate phase field simulations. In any case, the available computational 
resources will often set a limit on the smallest interfacial parameter s that one can compute 
for. Hence another aspect that needs to be taken into account in an objective comparison 
between phase field simulations and sharp interface approximations is the overall CPU time 
that is needed to obtain the results. While it can often be formally shown that phase field 
computations can attain an arbitrarily high accuracy, the existing limitations on computer 
hardware often mean that in practice very fine computations cannot be performed. In ad- 
dition, as discussed in [52j, the early computational approaches were limited as they could 
only be used in the presence of interfacial kinetics in the Gibbs-Thomson law. 

Historically these limitations of the phase field model have long been known, and as a 
result a different underlying interpretation of the model, the so-called "thin interface limit" , 
has been introduced and analyzed by Karma and Rappel [52i [53]. Their approach made it 
possible to do efficient computations with a smaller capillary length to interface thickness 
ratio, and to study the physically relevant case of small or zero kinetic coefficients. Later the 
findings of Karma and Rappel were reinterpreted as second order convergence with respect 
to the interfacial thickness, see |45l \32\ . 

The first successful phase field computations of dendritic growth were performed by 
Kobayashi [55j, and his computations demonstrated the importance of anisotropy for den- 
dritic growth. Since then many successful improvements with respect to numerical simu- 
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lations have appeared in the hterature. We refer to [HI |80l [52l [671 ES] and the references 
therein. 



FinaUy, we would hke to mention work on the numerical analysis of phase field and 
sharp interface approaches. Numerical analysis of discretizations of phase field models can 
be found in e.g. j29l SH EZl |33l EH [43l [I7j. Numerical analysis of discretizations of sharp 
interface models can be found in [THJ [12[ [13] . We also remark that level set methods are 
another possible way to solve Stefan problems and related free boundary problems. We 
refer to [65[ [72] for more details on how the level set method can be used to solve free 
boundary problems. 

The remainder of the paper is organized as follows. In Section [2] we state the sharp 
interface formulation of the two phase Stefan problem with kinetic undercooling and an 
anisotropic Gibbs-Thomson law. In Section [3] we state the corresponding phase field model 
and recall the finite element algorithms from [16j. In Section |4] we numerically compare the 
sharp interface method from [12j with the phase field algorithms from Section |3] for some 
isotropic benchmark problems with known true solutions. Computations for a phase field 
model with a correction term, for which a second order convergence property can formally 



be shown, are presented in Section 4.3, Finally, we compare the sharp interface and phase 



field methods for a variety of anisotropic model problems in Section [5) 



2 Sharp interface problem 

In this paper we concentrate on interfacial problems in materials science in which one driving 
force is due to capillary eflFects. In applications the interface often separates a solid and a 
liquid phase, say, or a solid phase and a gas phase. Let r{t) C M^, = 2,3, denote this 
sharp interface. Then the surface energy of r(t) is defined as 

7(n)ds, (2.1) 

m 

where n denotes the unit normal of T{t), and where the anisotropic density function 7 : 
]R>o with 7 G (7^(]R^ \ {0}) H C(]R^) is assumed to be absolutely homogeneous of 
degree one, i.e. 

7(Ap) = |A|7(p) VpGM^, VAgM ^ i{p) .p = -f{p) V p G \ {0}, 

with y denoting the gradient of 7. For all the considerations in this paper we assume that 
7 is of the class of anisotropics first introduced by the authors in ^ tUj; see also J2l Il6j . 



Relevant for our considerations is the first variation, — /^^, of (2.1), which can be com- 
puted as 

i^j := -V^ .7 (n) ; 

where V^. is the tangential divergence of F, see e.g. [30l [HI [12] . Note that /^^ reduces to 
the sum of the principal curvatures of F, in the isotropic case, i.e. when 7 satisfies 

l{p) = \P\ V J9 G (2.2) 
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2.1 Stefan problem 



Then the full Stefan problem we want to consider in this paper is given as follows, where 
r2 C M'' is a given fixed domain with boundary dfl and outer normal u. Find u : flx[0,T] ^ 
M and the interface {T{t))te[o,T] such that for all t G (0, T] the following conditions hold: 



du 



dn 



in 
-AV 



'dut-IC+Au = 



m 

/3(n) 
du 

m 



a 


To, 



au 



on On^^ u = ud 



in ^+{t) , 


(2.3a) 


on T{t) , 


(2.3b) 


on r{t) , 


(2.3c) 


on c^D^^ , 


(2.3d) 


in Q . 


(2.3e) 



In the above u denotes the deviation from the melting temperature Tm, i.e. Tm is the melting 
temperature for a planar interface. In addition, Q_{t) is the solid region, with boundary 
r(t) = dQ-{t)^ so that the liquid region is given by := ^\ Moreover, here 

and throughout this paper, for a quantity v defined on we use the shorthand notations 
V- := V \q_ and := v \q^. The parameters i9>0, A>0, p>0, Qf>0, a>0 
are assumed to be constant, while /C > is assumed to be constant in each phase. The 
mobility coefficient /3 : ^ ]R>o is assumed to satisfy I3(p) > for all p 7^ and to be 



positively homogeneous of degree one. We note that in the isotropic case (2.2) it is often 
also assumed that 

= \p\ V p G ^ /5(n) = 1 . (2.4) 

du-\- du- 
dn ~ dn 

T{t) in the direction of its normal n, which from now on we assume is pointing into Vtj^{t). 
Finally, dO. = dN^yJdo^ with dN^Hdo^ = 0, G M<o is the applied supercooling at the 
boundary, and Fq C and t^o ^ ^ ^ 1^ are given initial data. Here we use the convention 
that Ud = Q dO. = OnQ. 



In addition [/C |^]r(t)(^) •= (^+ — ^- ^^)(^) for all z G F(t), and V is the velocity of 



The model (2.3a-e) can be derived for example within the theory of rational thermody- 
namics and we refer to [48j for details. We remar k that a derivation from thermodynamics 
would lead to the identity ^ = t^- We note that (|2.3b) is the well-known Stefan condition, 
while (2.3c) is the Gibbs-Thomson condition, with kinetic undercooling if p > 0. The case 
i9>0,p>0,a>0 leads to the Stefan problem with the Gibbs-Thomson law and kinetic 
undercooling. In some models in the literature, see e.g. [58], the kinetic undercooling is set 
to zero, i.e. p = 0. Setting 1} = p = but keeping a > leads to the Mullins-Sekerka 
problem with the Gibbs-Thomson law, see [6IJ. 



We recall from [12j that for a solution u and F to (2.3a-e) it can be shown that the 
following equality holds 



^J'(F,ii) = -(/CVii,V^)-^ / 
dt a Jri 



d5 < 0^ 



(2.5) 
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where 



J^{T,u) := -\u-Ud\1 + — / -f{n)ds-XuD\n+{t)\ (2.6) 
^ ^ Jr(t) 

and where (•, •) denotes the L^^inner product over fi, with the corresponding norm given 
by I • |o, and where := /q^(^) 1 dx. 



2.2 Parametric method PFEM 

Traditional front tracking methods for sharp interface problems had a major drawback, in 
that the meshes used to describe the interface seriously deteriorated during the evolution. In 
addition, introducing mesh smoothing during the evolution is difficult, see e.g. the discussion 
in [To] . For interfaces in the plane it is possible to formulate a non-trivial method such that 
mesh points are nearly equally distributed during the evolution, see [50l [60] . The present 
authors introduced a novel parametric finite element method for problems involving curves 
and surfaces evolving in time, which has a simple variational structure and which leads 
to good mesh properties, see [8l [71 [TO]. In fact, for curves a semi-discrete variant leads to 
equally distributed mesh points in the isotropic case, while in the general anisotropic setting 
equidistribution with respect to some anisotropic weight function is obtained, see [9j. For 
surfaces the resulting meshes have also good properties, which in the isotropic case can be 
explained by using ideas from conformal geometry. In particular, no remeshing is needed 
during the evolution, even in the general anisotropic situation. An example triangulation 
obtained during the simulation of dendritic growth in three space dimensions can be seen 
in Figure [21} below. In addition, as the mesh for the parameterization of the interface is 
decoupled from the bulk mesh, no deformation of the bulk mesh is required in order to 
contain the interface at predefined locations on it. 



The novel and stable parametric finite element approximation of ( 2.3a| -e) in the case 



/C+ = JC- > has been introduced by the present authors in [12j, and this scheme has been 
extended to the more general case IC± > in [13j. Throughout this paper we will refer to 
these variants as PFEM. The algorithm PFEM features the discretization parameters /ir, 
hf^ he and r. Here hr refers to the fineness of the triangulated approximation of r(t), for 
which isoparametric piecewise linear finite elements are employed. In particular, a simple 
mesh refinement strategy allows for the natural growth of the interface, i.e. elements of the 
triangulated approximation of T{t) are refined when they become too large. Moreover, the 
temperature in the bulk is approximated with standard continuous piecewise linear finite 
elements, and hf and he refer to bulk mesh parameters for fine regions close to the interface 
and coarser regions far away from it. For all the computations presented in this paper we 
fix he = 8hf and, unless stated otherwise, we let hf ^ hr- Finally, r denotes a uniform 
time step size. The linear discrete systems of equations are solved with preconditioned 
conjugate gradient solvers of suitable Schur complement formulations. We refer to [T2[ [13] 
for more details. As indicated earlier, no remeshing of the discrete interface is necessary for 
the scheme PFEM, and all the numerical results presented in this paper for this scheme are 
performed without any redistribution of mesh points. 
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3 Phase field model 



We now state the phase field model that we are going to consider in this paper. To this 
end, for p G M^, let 



A{p) 



i{p) i{p) p 7^ , 

P = 0; 



and define 



\\x p = 0, 

where jl G ]R>o is a constant satisfying minp^o < A* < maxp^o 



(3.1) 



Moreover, let (p \ Q x (0,T) ^ M be the phase field variable, so that the sets {x G Q : 
±Lp{x^t) > 0} are approximations to with the zero level set of (f{-^t) approximating 

the interface T{t). On introducing the small interfacial parameter 6: > 0, it can be shown 
that ^ 

— Se{p) ^ / 7(n) ds, 

for e sufficiently small, where 

£e{^) := [ + £-^^{^) dx with c^:= [ ^2^{s) ds . (3.2) 

Jn J-i 

Here : M ^ [0, oc] is a double well potential, which we assume to be symmetric and to 
have its global minima at ±1. The canonical example is 



:= I (.2 _ 1)2 ^ v^'(,) 
Another possibility is to choose 

1^1 > 1 , 



and = ^ 2 2 



2 ' 



(3.3) 



(3.4) 



see e.g. [22l [23[ El ^40 J . Clearly the obstacle potential (3.4), which in contrast to the smooth 
quartic potential (3.3) forces (f to stay within the interval [—1, 1], is not differentiable at ±1. 
Hence, whenever we write ^^'(s) in the case (3.4) in this paper, we mean that the expression 
holds only for |5| < 1, and that in general a variational inequality needs to be employed. 



Our phase field model for (2.3a-e) is then given by the coupled system 
^wt + Xg{^)^t -V .{b{^)Vw) = rnQT •=^x (0,r) , 



w 



ud on do^ X (0, T) , 
on Dn^ X (0,r) , 
'dw{', Q) = dw^ in , 



(3.5a) 
(3.5b) 

(3.5c) 
(3.5d) 



7 



with 



a a 
dip 

du 



p 



in ill 



(3.6a) 







on dfl X (0, T) , (3.6b) 
in , (3.6c) 



where 



b{s) = |(l + s)/C+ + |(l-s)/C_ 
and where the function g G C-^(]R) is such that 



^(s)>0 Vse[-l,l] 



Q{y) dy = 1 and P(s) 



i?(?/)d|/. (3.7) 



We note that P, which is a monotonicaUy increasing function over the interval [—1,1] with 
P(— 1) = and P(l) = 1, is often caUed the interpolation function. In this paper, we follow 
the convention from where ^ = PMs called the shape function. Possible choices of g 
that will be considered in this paper are 



(i) q{s) 



ii) gis) = \ {l-s) , (iii) g{s) = If {s'-lf , (iv) g{s) = | (1-.^) . (3. 



More details on interpolation functions P, respectively shape functions g, can be found in 
e.g. [ini SSI EHj • In particular, if one also assumes that g is symmetric, i.e. 



gls 



g{-s) 



(3.9a) 



and that 

£)(!) = £)(-!) = 0, (3.9b) 
then a faster rate of convergence of the phase field model to the sharp interface limit, as 



6: ^ 0, can be shown in the isotropic case (2.2), (2.4) on replacing p in (3.6a) with the first 
order correction 

p:= p + epi, (3.10) 



where pi is defined in (4.7) in Section 4.3, below. The condition (3.9b) is one motivation 



for the latter two choices in (3.8), with the choice (3.8)(iv) also satisfying the stronger 



condition (4.9), below, for the quartic potential ( |3.3[). A n error analysis for a fully discrete 
approximation of the phase field model (3.5a-d), ( |3.6a|-c ) with the quartic potential (3.3) 
and the shape function (3.8)(i) in the isotropic case ( |2.2| ), ( |2.4D with OnQ = dQ and /C+ = 
JC- > has been performed in [43j. These authors also show convergence of the phase 
field discretizations to the underlying sharp interface problem as 6:, /i, r ^ 0, where h 
and r denote the discretization parameters in space and time, respectively. However, to 
our knowledge, no convergence rates are known for the convergence of discretizations of 
the phase field model to the sharp interface problem (2.3a e). Here we recall that for the 



much simpler situation of planar curvature fiow, as the sharp interface limit of the isotropic 
AUen-Cahn equation, such convergence rates have been obtained in [64j. In particular, it 
can be shown that the zero level sets of discretizations of 
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for the obstacle potential (3.4) converge with 0(e) to the sharp interface limit moving by 
V = /^if 

r = 0{h^) = 0{e^) . (3.11) 



While no such result is known for the full phase field model (3.5a-d), (3.6a c) even in the 
isotropic setting, it is natural to expect constraints of the form (3.11) in order to observe 
0(e) in practice. 



We remark that the phase field analogue of the sharp interface energy identity (2.5) is 
given by the formal energy bound 



w) = -(%) Vw,Vw)-e^ — (/x(V (^), ((^0') < 
at a 



(3.12) 



for the phase field model ( 3.5a| -d), ( |3.6a| -c), where 



_ . . . ,2 1 
:= - 1^ - Ud\o H 

Z Qj 



E,{if)-\UD I P((/^) dx. (3.13) 



Phase field models that satisfy such an inequality, in analogy to the sharp interface energy 
identity (2.5), are often called thermodynamically consistent, see [661 1791 US] . 



3.1 Phase field methods PF^^'-FEM and PF^^^-FEM 



Unconditionally stable, fully practical finite element approximations of ( 3.5a| -d), 
( 3.6a| -c) with either (3.3) or (3.4) have been introduced by the authors in p^. Here stable 
means that they satisfy a discrete analogue of the formal energy bound (3.12). Throughout 



this paper we will refer to the approximations from [16] for (3.3) and (3.4) as PF'^^^-FEM 
and PF^^^-FEM, respectively, where the inclusion of a subscript refers to the choice of shape 
function in ( |3.8 ), e.g. PF°^^-FEM. We recall from [16j that a side effect of the interpolation 
function P in (3.13) is that the function 



G{s) = a{ac^e) ^ "^{s) - ud^{s) 

need no longer have local minima at5 = ±lif2ix)7^0. This can result, for example, in 
undesired, artificial boundary layers for strong supercoolings, i.e. when —Ud is large. For 



the smooth potential ^ from (3.3), sufficient conditions for s = ±1 to be local minimum 
points of G{s) are 

q{±1) = q'{±1) = Q, (3.14) 



which is evidently satisfied by (3.8)(iii). In fact, in applications phase field models for 



solidification almost exclusively use this shape function; see e.g. [211 [3ll [59] . Pot the obstacle 



potential (3.4) the situation is similar, although there is more fiexibility in the possible 



choices of g. In particular, here a sufficient condition for G{s) to have local minima at 
5 = ±1 is given by 

a{ac^e)-^ ±udq{±1) > 0. (3.15) 
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On recalling that Ud < we see that for ( 3.15[) t o hold it is sufficient to require that 
^(1) = 0, which is evidently satisfied by (3.8)(ii), ( |3.8 )(iii) and (3.8)(iv). A major advantage 
of (3.8)(ii) over (3.8)(iii) and (3.8)(iv) is that for ( |3.8D (ii) it is possible to derive almost 
linear finite element approximations that are unconditionally stable. The corresponding 
unconditionally stable schemes for the nonlinear shape functions (3.8)(iii) and ( |3.8D (iv), on 
the other hand, turn out to be highly nonlinear. See [16j for more details. 



We remark that even when (3.14) and (3.15) are satisfied for the potentials (3.3) and 



(3.4), respectively, it is possible that mushy interfacial regions are observed in practice for 



approximations of the phase field model (3.5a-d), ( |3.6a| -c); see e.g. Figure 13, below. That 
is particularly the case in situations where the instability of the moving free boundary is 
strong, i.e. when —Ud aa~^ is large, recall (2.6) and see e.g. [61j. Then e needs to be chosen 



small, recall (3.13), so that the phase field variable ip admits well-defined interfacial regions 



that approximate the sharp interface T{t). This gives rise to a formal constraint of the form 



£<Ca{—UDCi) ^ if 2iD<0, 



(3.16) 



for the choice of the interfacial parameter e in terms of the physical parameters for the 



sharp interface problem (2.3a-e), irrespective of the choice of g. The reason for this is that 



in the estimate (3.12) the double well term ^(^) in Ss{(f) is for large s not strong 

enough to bound the unstable term involving which encourages the growth of the 

diffuse interface. 

The two algorithms PF^^'-FEM and PF^^^-FEM, which use continuous piecewise linear 
finite elements in space, feature the discretization parameters /i/, /ic and r. Here hf and he 
are mesh parameters for fine triangulations inside the diffuse interfacial region and coarser 
triangulations far away from it. Meaningful phase field simulations need to resolve the 
interfacial regions, whose width is of the order 5, and so a constraint of the form 



hf < Ce 



(3.17) 



needs to be enforced. For (3.4) the asymptotic interface width is tt^ in the isotropic case 
(2.2), and in this paper we always choose hf < ^ with hc = 8hf<7T£. Unless otherwise 
stated we let hf = Finally, r denotes a uniform time step size. Here we recall that 
the schemes PF^^'-FEM and PF^^^-FEM employ a semi-implicit discretization in time, 
which utilizes convex/concave splittings of the nonlinearity arising from the potential 
and from the interpolation function P. Such a splitting for was first proposed in [42j, see 
also PI, and the idea generalizes naturally to P; see Section 3.2, below, for details. This 



means that for the shape function choices (3.8)(i) and (3.8)(ii) almost linear schemes are 
obtained, while the choices (3.8)(iii) and ( 3.8| )(iv) give rise to more nonlinear finite element 
approximations; see [l6] for details. The discrete systems of linear equations and variational 
inequalities arising from the schemes PF^^^-FEM are solved with the Uzawa-multigrid solver 
from [4J , while the systems of nonlinear equations arising from PF^^^-FEM are solved with 
a Newton method. We refer to [16j for more details. 



For completeness we briefiy describe the choice of the initial profile (fo in (3.6c) in our 
numerical computations. Given the initial interface Fq from (|2.3e|), we let rfo ^ ^ ^ 1^ 
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denote the signed distance function of Fq. Then, on recaUing the asymptotic phase field 
profiles from e.g. [40 J, we define 

(po{x) = ^{e~'^do{x)) , where := { sin(s) |s| < | , (3.18a) 

5>f , 




for the obstacle potential (3.4), while for the smooth quartic potential (3.3) we use 

(fo{x) = ^{£~^ do{x)) , where $(5) := tanh(2"^ s) . (3.18b) 

For simplicity we use the profiles ( 3.18a| b) also in the anisotropic setting, where it would 
be more appropriate to replace do with a suitably defined anisotropic distance function rf^, 
see [36j for details. Finally, if i9 > 0, we fix wq = Uq. 



3.2 Possible time discretizations 



In Section |4] we will investigate the accuracy and the efficiency of several discretizations of 
the phase field model (3.5a-d), (3.6a^c) in the isotropic case (2.2), (2.4). In addition to the 
schemes PF^^'-FEM and PF^^^-FEM from fT6], which use a semi-implicit discretization in 
time, we will also look at a more implicit discretization and at a fully explicit discretization. 
For later reference, we now state the three different time discretizations, and for simplicity 
we do so on ignoring spatial discretization. A strong formulation of the time discretization 
from fl6j is given as follows. Let = + with being convex on M and 
being concave, and let P = P^ + P~ be a similar splitting that is convex/ concave on a 
suitable superset of [—1, 1], where we recall that (f need not remain in [—1, 1] for the quartic 
potential (Q. We also define := (P^)^ For the schemes PF^^^-FEM and PF^^^-FEM 



we set (5) = — I 5^ and 
(i) g+{s)=0, 



(ii) g+{s) = 0, (iii) g+{s) 



2 ^ ' 



2s 



for the choices of ^ in (3.8). The semi-imphcit time discretization employed by the schemes 



Ppobs.pEM and PF-^^^-FEM can then be formulated as: 

(Wn - Wn-l) + A (g'^iiPn) + g~ {(fn-l)) {^n " ^n-l) " T V . {h{ipn-l)^ Wn) = , (3.19a) 

e ^ ^" ~ ^"-^ -e^^n + e-' ((*+)'((^n) + (*-)'(¥'n-i)) . 

a a r 

(3.19b) 



With this time discretization existence of a unique solution [lp^^w^) to the fully discrete 
scheme can be shown for arbitrary time step sizes r if ^+ = 0, where may not be 
unique if i9 = and = in very rare circumstances. Moreover, any solution to the 
semi-implicit schemes PF^^^-FEM and PF*^^^-FEM is stable; see [16j for details. The semi- 
implicit time discretization in ( 3.19a| b) can be modified to an implicit time discretization 
by replacing ( |3.19bD with 



a 



£ 

a T 



(3.20) 
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which then gives rise to a time step size constraint of the form 



T < pa ^ 



if p > , 



(3.21) 



in order to ensure the existence of a unique, stable solution in the case ^+ = and p > 0. 
In the situation i9 = p = and q{^s) = |, with = dQ^ which has been treated by the 
present authors in [15j, a stronger time step constraint of the form r = 0{£^) arises for 
the implicit discretization (3.20); see also [221, '20j. Here we recall that it is often observed 



that implicit time discretizations of AUen-Cahn and Cahn-Hilliard type equations yield a 
better accuracy in time compared to semi-implicit time discretizations as in ( |3l9al b); see 
e.g. [20l [211 HTj. We will present several computations for a variant of PF^^^-FEM with the 



implicit time discretization (3.19a), (3.20) in Sections and ^ 



Finally, fully explicit approximations, as advanced in e.g. 



sidered. Here we replace (3.19a,b) by 



n-1 



n-1) 



0, 



can also be con- 



(3.22a) 



a , 
a 



-l)^n-l = ^ 



a r 



- e A + + (^n-i)) . (3.22b) 



If i9 > and p > 0, then the above fully explicit time discretization is well-defined, and in 
this case stability of the fully discrete scheme can be shown if 

r = 0{h^) . (3.23) 

In particular, in the case of the obstacle potential (|3.4|), if r < \'d e p{\^ a)~^ g'^^^^ then 



the solutions to the fully discrete variant of ( |3.22a ,b) are stable if 

r/i"^ < min{| i9/Cj^^x, pa"-^} 



(3.24) 



where ^^ax niax5^[_i^i] 1^(^)1, /Cmax max{/C+,/C_} and where is a constant only 
depending on the spatial mesh. The advantage of ( 3.22a| b) over ( 3.19a| b) is that the 
discretized systems of equations now decouple in space, which leads to huge efficiency gains 
when the computations are performed in parallel on a large cluster. However, in practice this 



advantage is often negated because (3.23), together with e.g. (3.17), enforces that very small 
time steps need to be taken. In Section |4.2| and in Section [5 we will present computations 
for a variant of PF^^^-FEM with an explicit time discretization as in (3.22a,b). 



4 Quantitative comparison for isotropic problems 

A standard validation used for phase field models in the literature is the comparison of tip 
velocities between the computed phase field discretizations and real world measurements 
from the laboratory, see e.g. ^2l|67j. Here the physical parameters in the phase field model 
have to be chosen appropriately, so that they correspond to the physical properties of the 
material in question. However, often the exact values of these parameters are unknown or 
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Figure 1: The true solution from (4.2a) at times t = 0, 10 ^, 2 x 10 ^. 



they are themselves based on measurements. Here we propose a much simpler quantitative 
validation, which makes use of known radially symmetric solutions to the underlying sharp 
interface problem in the isotropic case. We would argue that such a simple comparison 
should be part of the validation of every phase field method to be proposed in the literature. 
It gives an indication of the accuracy of the overall method and it helps to fine-tune the 
discretization parameters that should be used for the anisotropic physical applications. 



In particular, in this section we consider the following isotropic variant of ( |2.3a| -e). Find 
: X [0, T] ^ M and the interface (r(t))tG[o,T] such that for all t G (0,T] it holds that: 



'dut — Au = f 
du 



-V 



m 

du 
dv 

m 



a — u 



= on On^^ u = ud 
= ro, ^u{-,0)=^Uo 



in Q \ r{t) , 


(4.1a) 


on r{t), 


(4.1b) 


on T{t), 


(4.1c) 


on do^ , 


(4.1d) 


in Q . 


(4.1e) 



Here / : [0, T] ^ M in (|4.1aD is a given spatially homogeneous forcing term. In the phase 
field approximation (3.5a-d), ( |3.6a| -c) this forcing appears analogously as a right hand side 
term / in (3.5a). Note that for / = the above system ( 4.1a| -e) corresponds to ( 2.3a| -e) 
with ([^21), pl| and /C± = a = A = 1. 



4.1 Mullins-Sekerka problem 

We first consider the quasi-static case = p = 0. To this end we take the known solution 
of an annular region f^-(t), for which the inner boundary shrinks to a point so that Q-(t) 
becomes a disk for sufficiently large t. Here we take a = 1 and let = dQ. In addition, 
r(0) = Fq consists of two concentric circles/spheres. It is then not difficult to show that 
the two radii ri < r2 satisfy the following system of nonlinear ODEs: In the case d = 2 we 
have 

N* = --i^ and [r2], = -li^ = ^[ri], VtG[0,To), (4.2a) 
ri In — To m — ro 

^ ri ^ ^-j^ z, 
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0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 ^"^ 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 
true solution k=5 k=6 k=7 true solution k=5 k=6 k=7 




1.5' ' ' ' ' ' ' ' ' 1.5' ' ' ' ' ' 1 1 ' 

0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 

true solution k=5 k=6 k=7 true solution k=5 k=6 k=7 

Figure 2: (PFJ^^-FEM, e'^ = Stt, 16 tt, 32 tt, 64 tt) Comparison of the energies T and 
for the benchmark problem [T] with T = 4 x 10~^. The uniform time step sizes are chosen 
as r = 10-^ = 5 ^ 7. 



while for = 3 it holds that 
2 ri + r2 



[ri]t = 



ri 



and [r2\t 



4^^ = In* vtG[o,ro), 

r| r2 — ri r| 



(4.2b) 



where Tq is the extinction time of the smaller sphere, i.e. limt^To ^i(^) = 0, see e.g. [HI [76]. 
Note that the corresponding solution u satisfying ( 2.3a| ^e) is given by the radially symmetric 
function 



u{x^ t) < 



d-i 



In 



r2(t)-ri(t) 



\x\ ri(t) ^ r2(t) 
ri(t) 
\x\ r2(t)-ri(t) 



d = 2 



d 



d-i 



\x\ > r2{t) 



ri{t) < \x\ < r2{t) 



(4.3) 



As (4.2a,b) does not appear to be analytically solvable, it needs to be integrated numerically 
to compute the solution (ri,r2)(t), for t G [0,T], where T < Tq. Possible strategies to 
integrate ( 4.2a| b) to a high accuracy are described in [l2]. We visualize the evolution of 
r(t) over the time interval [0, 2 x 10~^] in Figure [l] The above true solution forms the basis 
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0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 
true solution k=5 k=6 k=7 




0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 
true solution k=5 k=6 k=7 



0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 
true solution k=5 k=6 k=7 



Figure 3: (PF^^^-FEM with the imphcit time discretization from (3.20), s ^ = 
Stt, 16 tt, 32 tt, 64 tt) Comparison of the energies and J^^ for the benchmark problem [l] 
with T = 4 X 10~^. The uniform time step sizes are chosen as r = 10~^, k = 5 ^ 7. 



of our first benchmark problem: 



Benchmark problem 1: 2d Mullins-Sekerka with = p = 0. 



True solution (4.2a), (4.3) to (4.1a-e) with i9 = p = and a = 1. 
Initial data (ri,r2)(0) = (0.24,0.4). 
Domain = (— |, with dQ. 

Time interval [0,r] with T = 10'^ so that (ri,r2)(r) ^ (0.18,0.37). 



In Figure [2] we compare the energy of th e tru e sharp interface solution, recall (2.6), 
to the corresponding energies J^^ ^ T^^ recall (3.13), of the finite element approximations 
from the algorithm PF^^^-FEM on the time interval [0,4 x 10~^]. Here we recall that for 
the given data and the given evolution it holds that 



j'(r,^) = |r(t)| 



Id5 = 2 7r(ri(t)+r2(t)). 



Very similar energy plots can be obtained for the other variants of PF^^^-FEM and those of 
PF°^^-FEM. We note that for decreasing 5, the time step size r needs to be chosen smaller 
and smaller in order to capture the correct time scaling of the evolution. We compare this 
computation for PF^^^-FEM, which uses a semi-implicit discretization in time, now with 



one computation for the implicit time discretization as in (3.20). Here we recall that better 
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1.5 ' 

0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.0O4 

true solution i=4 i=6 

j=3 1=5 1=7 

Figure 4: (PFEM) Comparison of the energies T and for the benchmark problem [l] 
with T = 4 X 10~^. The uniform time step size is chosen as r = 10~^, while the spatial 
discretizations are proportional to 2~-^ , j = 3 ^ 7. 



accuracy for such discretizations has been reported in [20] for the isotropic Cahn-Hilliard 
equation, i.e. (3.5a), (3.6a) for (2.2) with (|3.4D, (3.8)(i) and with i9 = p = 0. In fact, here we 



observe a similar behaviour. See Figure [3| where even for very large choices of r the time 
evolution of the phase field model seems to be captured accurately. Finally, we show some 
discrete energies ^ T from the parametric scheme PFEM in Figure [i} Here we choose 
rather crude discretization parameters, since otherwise the discrete energies would lie 
virtually on top of the true energy T . 

In Table pi we present the errors in r\ and in u for the scheme PF^^^^-FEM for a selec- 
tion of interface parameters e and for a range of discretization parameters hj and r for 
the benchmark problem [!} The displayed error quantities are defined as 



Eo<n<T/r - ^(•.^^)lli2(^))^ and 

ri(nr)|, where r^^(t) := infjs > : 99^(5 ei,t) 



\w — u\\i2 := 
rillLoo := maxo<n<r/r - 
0}, with ei = (1,0)^ being the first unit 



vector in M , denotes the phase field approximation of the inner radius. We also show the 
number of degrees of freedom (DOFs) for the calculation of the discrete solution for the final 
time step at time t = T. The presented overall CPU times are for a single-thread computa- 
tion on an Intel 17-860 (2.8 GHz) processor with 8 GB of main memory. For the benchmark 
problem [1] the remaining variants of PF^^^-FEM and PF^^^-FEM exhibit very similar errors 
to the ones in Table [T} and so we do not present them here. In later computations we 
will also employ the stronger norm — 2i||7,oo := maxo<n<r/r \\w^{-^nT) — nr)||i,oo(^) 
for the temperature error. However, for the experiments in Table [T] no convergence can be 
observed in the L^(J77^)-error for the true temperature (4.3) for the phase field approxima- 
tions. In fact, for the computations in Table [T the errors \\w^ — u\\loo are in the interval 
[1, 16], where we note that the true solution (4.3) itself remains in the range [—2.8, 5.6] over 
the computed time interval. It is for this reason that we report the weaker error norms 
11'^^ ~ '^||l2 in Table [ij 

A repeat of the computation in Table [T} but now for an implicit time discretization of 
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£-1 


2^^/ 


T 


Ik^i - tiIIloo 


\\w'^ — u\\l2 


DOFs(T) 


CPU time 


Stt 


64 


10" 


-5 


1.3004e-02 


1.8448e-03 


4522 


5 sees 




64 


10" 


-6 


1.3154e-02 


9.9918e-04 


4402 


39 sees 




64 


10" 


-7 


1.3175e-02 


3.2861e-04 


4354 


6:54 mins 




128 


10" 


-7 


1.4283e-02 


3.5051e-04 


15626 


37:22 mins 


167r 


128 


10" 


-5 


1.4139e-02 


2.1678e-03 


10378 


15 sees 




128 


10" 


-6 


2.6476e-03 


2.5879e-04 


10130 


1:58 mins 




128 


10" 


-7 


4.9843e-03 


1.0553e-04 


10106 


17:06 mins 




256 


10" 


-7 


5.5154e-03 


1.1024e-04 


34682 


1:23 hours 


32 TT 


256 


10" 


-5 


3.4984e-02 


5.6999e-03 


23362 


50 sees 




256 


10" 


-6 


6.2022e-03 


3.4521e-04 


21650 


6:28 mins 




256 


10" 


-7 


8.8958e-04 


3.7115e-05 


21082 


50:03 mins 




512 


10" 


-7 


1.8543e-03 


3.3740e-05 


74322 


5:04 hours 


64 TT 


512 


10" 


-5 


5.0345e-02 


8.2731e-03 


52602 


1:54 mins 




512 


10" 


-6 


2.2346e-02 


1.2118e-03 


49370 


17:19 mins 




512 


10" 


-7 


4.7000e-03 


1.1020e-04 


46922 


2:14 hours 




1024 


10" 


-7 


2.1723e-03 


3.5938e-05 


166498 


15:50 hours 


128 TT 


1024 


10" 


-5 


5.6324e-02 


9.1597e-03 


122314 


7:25 mins 




1024 


10" 


-6 


4.2515e-02 


2.2390e-03 


118818 


1:01 hours 




1024 


10" 


-7 


1.5232e-02 


2.8138e-04 


112794 


9:50 hours 




2048 


10" 


-7 


1.0919e-02 


1.8112e-04 


404962 


68:23 hours 



Table 1: Benchmark problem fTl for PFg'-FEM. 



PF?.r-FEM can be seen in Table 



(i) 



and 



One clearly observes that, for fixed 6:, the errors ||r: 



XI 



\W 



u\\l2 soon appear to be almost independent of the time step size r. This 



indicates that the implicit time discretization from (3.20) manages to eliminate the temporal 
discretization error relatively quicker than the semi-implicit discretization from (3.19a,b). 
Moreover, for small £ and fixed r, the error in the approximation of the sharp interface 
problem ( 4.1a| -e) is in general significantly smaller for the implicit time discretization. We 
remark that the converged errors in Table [2| appear to indicate a convergence of 0{£) in the 



error ||r^^ — ri||i,oo, with a similar convergence rate for the error H'u;^ — 'uHl^, if discretization 
errors are neglected. 

Finally we note that the numbers in Tables [T] and [2] indicate that refining the mesh 
discretization parameters /i/, and hence /ic, in general does not reduce the error. Hence 
choosing he = 8hf = n s appears to be sufficient for classical phase field model computa- 
tions, and we will restrict ourselves to this choice from now on in this paper. 

We compare these convergence experiments with the corresponding errors for the sharp 
interface algorithm PFEM in Table |3j For these sets of experiments we always choose 
hr ^ hf. Here the error quantities \\u^ — u\\loo and \\u^ — u\\l2 are defined as \\w^ — u\\loo 
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e-1 


2V/1/ 


r 


Ik^i - tiIIloo 


\\w^ — u\\i,2 


DOFs(T) 


CPU time 


Stt 


64 


10" 


-5 


3.5385e-02 


7.8725e-03 


5018 


10 sees 




64 


10" 


-6 


3.5228e-02 


2.4832e-03 


5018 


1:06 mins 




64 


10" 


-7 


3.5211e-02 


7.8429e-04 


5018 


12:30 mins 




128 


10" 


-7 


3.6235e-02 


7.9552e-04 


18658 


50:07 mins 


167r 


128 


10" 


-5 


1.4453e-02 


2.2343e-03 


11514 


22 sees 




128 


10" 


-6 


1.4667e-02 


7.2292e-04 


11530 


2:20 mins 




128 


10" 


-7 


1.4675e-02 


2.2894e-04 


11514 


27:29 mins 




256 


10" 


-7 


1.5018e-02 


2.3527e-04 


38858 


1:31 hours 


32 TT 


256 


10" 


-5 


5.6651e-03 


8.2274e-04 


24802 


58 sees 




256 


10" 


-6 


6.4112e-03 


2.9604e-04 


24290 


7:23 mins 




256 


10" 


-7 


6.4838e-03 


9.6155e-05 


24322 


1:06 hours 




512 


10" 


-7 


6.6777e-03 


1.0156e-04 


83690 


5:41 hours 


64 TT 


512 


10" 


-5 


3.1282e-04 


1.3540e-04 


54298 


2:20 mins 




512 


10" 


-6 


2.5224e-03 


1.1540e-04 


53466 


17:31 mins 




512 


10" 


-7 


2.8135e-03 


4.4061e-05 


53114 


3:27 hours 




1024 


10" 


-7 


3.0098e-03 


4.2577e-05 


185410 


20:03 hours 


128 TT 


1024 


10" 


-5 


8.3850e-03 


1.4662e-03 


128866 


15:56 mins 




1024 


10" 


-6 


3.3622e-04 


2.2218e-05 


123770 


1:17 hours 




1024 


10" 


-7 


8.3440e-04 


1.5965e-05 


123402 


10:46 hours 




2048 


10" 


-7 


1.1768e-03 


1.6903e-05 


437330 


116:48 hours 



Table 2: Benchmark problem jlj for PF^^^-FEM with the implicit time discretization from 



and \\w^ — u\\l2 as before, but with replaced by v!^. In addition, we let \\r\ — ri||i,oo := 
maxo<ri<r/r ^^^pGrf(nr) lli^l "^il^"^)!? with T\{t) denoting the parametric approximation of 
the inner circle of the true solution T{t). Note that the norm in the definition of ||r5^ — ri H^oo 
employed here is much stronger than the phase field equivalent ||r^^ — ri||i,oo introduced 
earlier, where the difference between the true interface position ri{t) and the phase field 
approximation is measured in the Xi-coordinate direction only. All of the error quantities 
shown in Table [3] appear to be converging with order at least 0{h) if the time discretization 
errors are neglected. 

[3] convey a very clear message. Firstly, we recall that the 



2| do not converge in the norm \\w^ — u\\loo^ whereas \\u^ 



The numbers in Tables [T 
experiments in Tables [l] and 

1^11^00 in Table [3] does appear to converge with 0{h). Secondly, we can see that even with 
computations that take almost 5 days, the phase field schemes PF°^^^-FEM and PF^^^-FEM 
cannot reduce the error in the radius to below 3 x 10""^. Yet, better accuracies for the 
radius can be achieved with the sharp interface approximation PFEM, running on the same 
computing hardware, in less than a minute. Hence, for this measurement, the simulations 
with PFEM are at least 7000 times faster than the computations with PF^^^-FEM and 
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T 


Iki - nlU"" 


\\u'^ — u\\loo 


— u\\l2 


DOFs(r) 


CPU time 


8 


10" 


-4 


3.6782e-02 


2.8260e-00 


3.5916e-02 


113 


sees 


8 


10" 


-5 


4.2433e-02 


4.8862e-00 


4.1925e-02 


113 


sees 


16 


10- 


-4 


1.1028e-02 


9.5419e-01 


1.2123e-02 


285 


sees 


16 


10" 


-5 


1.3394e-02 


1.7413e-00 


1.5531e-02 


285 


sees 


32 


10" 


-4 


3.6288e-03 


4.3183e-01 


4.8601e-03 


693 


sees 


32 


10" 


-5 


5.7298e-03 


6.3045e-01 


7.1794e-03 


657 


1 sees 


64 


10" 


-4 


7.7318e-04 


2.9301e-01 


2.8734e-03 


1585 


sees 


64 


10" 


-5 


2.4803e-03 


2.9384e-01 


3.2140e-03 


1473 


1 sees 


128 


10" 


-4 


5.8266e-04 


3.1247e-01 


3.2813e-03 


3553 


sees 


128 


10" 


-5 


1.1406e-03 


1.1283e-01 


1.4618e-03 


3213 


3 sees 


128 


10" 


-6 


1.3262e-03 


1.5152e-01 


1.7016e-03 


3173 


32 sees 


256 


10" 


-4 


1.1384e-03 


3.4935e-01 


3.6955e-03 


8289 


1 sees 


256 


10" 


-5 


4.7583e-04 


5.9495e-02 


6.7789e-04 


6945 


7 sees 


256 


10" 


-6 


6.4981e-04 


7.6398e-02 


8.4454e-04 


6777 


1:11 mins 


512 


10" 


-4 


1.4098e-03 


3.6302e-01 


3.8770e-03 


20381 


3 sees 


512 


10" 


-5 


1.3003e-04 


3.4588e-02 


3.4619e-04 


16109 


19 sees 


512 


10" 


-6 


2.9886e-04 


3.4579e-02 


3.9724e-04 


15649 


3:02 mins 


1024 


10" 


-4 


1.4724e-03 


3.6426e-01 


3.9180e-03 


43133 


6 sees 


1024 


10" 


-5 


3.8716e-05 


3.1143e-02 


3.2082e-04 


41757 


57 sees 


1024 


10" 


-6 


1.3093e-04 


1.4765e-02 


1.8090e-04 


39493 


9:31 mins 


2048 


10" 


-4 


1.5121e-03 


3.6991e-01 


3.9402e-03 


93385 


18 sees 


2048 


10" 


-5 


1.1372e-04 


3.7037e-02 


3.6346e-04 


120429 


3:35 mins 


2048 


10" 


-6 


5.3291e-05 


6.9988e-03 


7.7964e-05 


112285 


32:22 mins 



Table 3: Benchmark problem [T] for PFEM. 



PF^^^-FEM. The main reason behind this very slow convergence appears to be that the 
biggest contribution to the observed error comes from the interfacial parameter s. Hence 
in order to obtain reasonable errors, s needs to be taken very small, which on recalling 



(3.17) implies that the discretization parameters need to be chosen very small as well; recall 



also (3.11). Unfortunately, phase field computations thus soon reach the limit of what is 



computable on today's computer hardware. As an aside we note that when comparing CPU 
times between e.g. PF^J^^-FEM and PFEM in terms of degrees of freedom, then it is crucial 
to take into account the value of r, as this will be indirectly proportional to the number 
of algebraic systems that need to be solved during the whole simulation. In fact, for the 
same number of degrees of freedom and the same value of r, the CPU times between the 
two different algorithms are similar. 

To better visuahze the relative performances of PF^J'-FEM in Table |l[ PFj^^'-FEM in 
Table [2] and PFEM in Table [3| we present a plot of the errors in the radms ri against the 
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Figure 5: (Benchmark problem [I]) Log-log scatter plot of ||r^^ — riH^oo and \\r^ — riH^oo 
against the CPU time for the entries in Table fTl (PFg'-FEM, blue rhombi), Table [2| (PFj^^- 
FEM, green rhombi) and Tablels] (PFEM, red circles). 



necessary CPU time for all the entries in the three tables in Figure [5) This plot underlines 
the superiority of the sharp interface algorithm PFEM over the phase field methods. 



For the 3d solution (4.2b), (4.3) it is virtually impossible to perform a meaningful 
convergence experiment for the phase field approximations PF^^^-FEM and PF^^^-FEM. 
The reason is that for the same values of e and for comparable discretization parameters 
as in e.g. Table [T} the simulations in three space dimensions do not finish in a reasonable 
amount of time. For an example of a convergence test for the solution (4.2b), (4.3) for the 
scheme PFEM we refer to ^ Table 6]. 



4.2 Stefan problem 



In this section we consider the full Stefan problem (4.1a-e) with = 1 and p > 0. Here we 



adapt the expanding circle/sphere solution introduced in [70^ p. 303-304], where the radius 
of the circle/sphere is given by r{t) with 



r{t) = {r^{0)+ty 



In particular, we let 

z{t) = 

Then it is easy to see that 



a{d — 1) + I p 



v{s) = 



e 4 



1 y 



- dy. 



u(x, t) = 



z{t) 



\x\ < r{t) , 



(4.4a) 



(4.4b) 
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is the solution to (4.1a-e) with 

m 



d 



z{t) 



a{d — 1) 



and with ud in (4. Id) replaced by u\d^Q on = dVt. 



(4.5) 



For the Stefan problem with kinetic undercooling we propose the following benchmark 



problem, where on recalling e.g. (3.16) we note that increasing the parameter ^ G N leads 



to the benchmark problem becoming computationally more and more challenging. 
Benchmark problem 2^^^: 2d Stefan problem with i9 = 1 and p > 0. 



and with d = 1 and a = p = 10 



True solution (4.4a,b) to (4.1a-e) with (4.5 
Initial data r(0) = \ and = u{-^ 0) from (4.4b) 
Domain = (—1, 1) with = and = u\o^q from (4.4b). 
Time interval [0,r] with T = |, so that r{T) = |. 



For the benchmark problem 2^^^ with ^ = 1, all of the phase field schemes were able 
to compute the necessary evolutions reasonably well. As an example we show the results 

is the same as 



for the scheme PF^^^-FEM in Table |4| Here the definition of ||r^^ 
W^xi ~ ^i|U°° with ri replaced by r. In order to be able to assess the absolut e tem perature 
errors \\w^ — u\\loo^ we note that for this benchmark problem the true solution (4.4b) remains 



in the range [—0.95, —0.15]. The same computation now for the implicit time discretization 
from (3.20) is shown in Table [s] In order to visualize the different behaviour of the two 
diflFerent time discretizations, we plot the scaled phase field approximations £^ ^ recall 
(3.2), together with |r(t)| = 27rr{t) in Figures |6] and [?} Similarly to Section 4.1 it can 
be seen that the implicit time discretization eliminates the time discretization error in 
the approximation of the phase field equations sooner than the semi-implicit discretization. 
Moreover, for the smallest value of s the absolute errors ||r^^ — r H^oo appear to be significantly 
smaller for the implicit scheme. We remark that the converged errors in Table [5] appear to 
indicate a convergence of 0{£) in the error ||r^^ — with a similar convergence rate for 

the error — t^||i,2, if discretization errors are neglected. 



The results for the same benchmark problem for the scheme PF^^J-FEM with the explicit 
time discretization from ( 3.22a| b) are shown in Table [6j Here the reported CPU times 
are given as guidelines only, because we did not employ any parallelization. Recall that 
the discrete systems of equations decouple in space, and so a significant speedup of the 
computations can be expected if they are run in parallel on a large cluster. Firstly, we see 
that the obtained results appear to confirm the stability constraint ( |3.24 ), i.e. r <\C^h^ . 
Secondly, it can be observed that once the explicit method is stable, there is hardly any 
variation in the numerical results when decreasing r further. And finally, it is clear from 
Table [6] that there is no convergence in the reported error quantities, i.e. the phase field 
simulations do not converge to the sharp interface problem (4.1a-e) as 6:, /i/, r 0. We 
conjecture that this phenomenon is due to the sensitivity of the explicit method to the 



employed mesh adaptation strategy, recall Section |3.1[ This is confirmed by repeating the 
simulations for the explicit scheme on uniform grids, see Table [7| Now the errors appear to 
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Table 4: Benchmark problem 



2W 



with ^ = 1 for PF^.^.^-FEM. 

(ill) 



be converging, and the absolute errors agree with the corresponding converged errors from 
the semi-implicit and implicit variants of PF^??J-FEM, on which we do not report here. 

We compare the above phase field errors with the corresponding errors for the sharp 
interface algorithm PFEM. Here ||r^ — rH^oo is defined as HrJ^ — riH^oo, but with F^ replaced 
by F^ and with ri replaced by r. The errors for the benchmark problem 2^^^ with i = 1 
are reported in Table |8| Comparing the results in Tables [4]-[8] reveals once again that 
the sharp interface approximations from PFEM are more accurate than the corresponding 
computations from the phase field schemes PF^^^-FEM and PF^^^-FEM, and they can be 
obtained in a fraction of the CPU time. For example, we see from the Tables |4| [5] and [8] that 
in order to reduce the error in both F and u below 10~^ requires 13 seconds with PFEM, 
but it takes around 5 and 11 hours, depending on the time discretization, with PF^^^-FEM. 
In other words, in this measure the parametric front tracking method PFEM is between 
1 300 and 2 900 times faster than the phase field methods. A visualization of the numerical 



results in Tables HHS^ can be seen in Figure 10, below. 



If we increase the parameter i in the benchmark problem 2^^^ to ^ = 3, which on 



recalling (3.16) means that the problem is now computationally more challenging, then 
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cretization from (3.20) 



Table 5: Benchmark problem 2^^^ with ^ = 1 for PF^^^-FEM with the implicit time dis- 



for moderate values of e all of the phase field schemes exhibit mushy regions in which the 
phase field approximations Lp^ in modulus take on values significantly smaller than unity. 
This leads to excessive CPU times, since the adaptive mesh strategy uses fine meshes in 
the interfacial regions. In addition, often the mushy interfacial region quickly reaches the 
external boundary dVt^ which creates additional interfaces, so that the phase field solutions 
Lp^ no longer approximate the radially symmetric sharp interface solution V{t). Hence in 
what follows we present convergence experiments for the benchmark problem 2^^^ with ^ = 3 
only for > 16 tt. See Table jo] for the results for the scheme PF^j*?J-FEM, where we note 
that for this benchmark problem the true solution (4.4b) remains in the range [—0.35,0]. 
Similar results can be obtained for the other variants of the schemes PF^^^-FEM and PF^^^- 
FEM that satisfy ^(1) = 0, but we omit them here for brevity. In particular, the scheme 
PF?J??-FEM with the implicit time discretization from (3.20) yields almost identical results 



(ii)" 



to the ones in Table |9| for the step sizes r=10 ,/c = 5^6. We remark that for the 
numbers in Table [9] it is somewhat speculative to infer convergence rates in terms of 6:, since 



the errors 



XI 



and 



converged yet in terms of h 



\w — u 
and r. 



I LOO for the smallest value of £ do not appear to have 
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Figure 6: (PFJ^-FEM, e'^ = 47r, Stt, 16 tt, 32 tt) Comparison of |r(t)| and c'^ for the 
benchmark problem 2^^-^ with i = 1. The uniform time step sizes are chosen as r = 10~^, 



We again compare these numbers with the corresponding errors for the sharp interface 



algorithm PFEM, see Table [10} The results in Tables [Q] and [TOl confirm once more that the 
sharp interface approximations from PFEM are more accurate. For example, in order to 

M, 



reduce both the error in F and in u to below 10~^ requires about a minute with PF 
but 46 minutes with PF^??J-FEM. But crucially, it is clear from the numbers in Tables 
and [9] that only by decreasing s further can the observed errors for the phase field metKods 
be reduced. This in turn will lead to enforced reductions in hf and r, recall e.g. (3.11) 



llEI 



(3.17), (3.21) and (3.23). Overall this makes it impossible to perform these computations in 



practice. On the other hand, the presented errors in Tables [8| and \T0\ for the scheme PFEM 
indicate a convergence in the error ||r^ — rH^oo of order at least 0{h)^ if time discretization 
eflFects are neglected. Apart from the run for the finest value of /i/ in Table [8} where the 
time discretization error does not seem to have been eliminated yet, the same can be said 
about the temperature error 



— u\ 



LOO . 



In order to visualize the relative performances of PF^??^^-FEM in Table |9| and PFEM in 



Table [T0| we present a plot of the errors in the radius r against the necessary CPU time for 
all the entries in the two tables in Figure [8| 



Similarly to Section |4.1[ it is not possible to perform a meaningful convergence test for 
the solution (4.4a,b) in the case = 3 for the phase field approximations PF^^^-FEM and 
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(PF^^^-FEM with the imphcit time discretization from (|3.20|), e' 



Figure 7: 

47r, Stt, 16 tt, 32 tt) Comparison of |r(t)| and c^^ for the benchmark problem 2^^^ with 
£ = 1. The uniform time step sizes are chosen as r = 10~^, k = 2 ^ 6. 
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Figure 8: (Benchmark problem 2^^^ 



with ^ = 3) Log-log scatter plot of ||r^^ — r||i,oo and 
against the CPU time for the entries in Table |9] (PF^^?^^-FEM, blue rhombi) and 
Table fTU](PFEM, red circles). 
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Table 6: Benchmark problem 
cretization from ( 3.22a| b). 



with ^ = 1 for PF^^)'-FEM with the explicit time dis- 



PF^^^-FEM. As an example for such a convergence experiment for the parametric scheme 



PFEM we now consider the three-dimensional analogue of the benchmark problem 2^^^ 



Benchmark problem 3^^^: 3d Stefan problem with d = 1 and p > 0. 


Same as benchmark problem 


2W 


, but on the domain = (— 1, 1)^. 



The corresponding errors are shown in Table |TT| where we let /ir 



6hf. Similarly to 

the two dimensional benchmark problem, the two errors ||r^ — r||i,oo and — t^||i,oo appear 
to converge with order at least 0{h)^ if time discretization effects are neglected. 

Our final benchmark problem is the Stefan problem without interfacial kinetics in the 
Gibbs-Thomson law. Here we recall from e.g. [52j that often standard, classical phase field 
methods are not able to deal with this case in practice. 



Benchmark problem 4^^^: 2d Stefan problem with = 1 and p = 0. 


Same as benchmark problem 


2W 


, but with a = 10 ^ and p = 0. 



We stress that the approach from [16j for the phase field system has no problems in 
dealing with the case without interfacial kinetics, i.e. if p = 0. For example, a computation 
for the scheme PF???^^-FEM can be found in Table [l2l We compare these results with the 



(ii) 



corresponding computation for the sharp interface algorithm PFEM in Table 13, where 
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Table 7: Benchmark problem 
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cretization from ( 3.22a| b) and wit 
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Fi gure 9: (Benchmark problem with ^ = 3) Log-log scatter plot of ||r^^ — r||i,oo and 



|r — r ll^oo against the CPU time for the entries in Table 



Table 13 (PFEM, red circles). 
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(PF^^?)'-FEM, blue rhombi) and 



once again it appears that the error quantities converge with 0{h) if the time discretization 
errors are neglected. We visualize the relative performances of PF^^J-FEM in Table 12 and 
PFEM in Table [13] in Figure [9} As before, the performance of PFEM is vastly superior to 
the corresponding phase field computations. 



4.3 Second order accurate isotropic phase field model 



In this subsection we recall a variant of the phase field model ( 3.5a| -d), (3.6a~c) which in 
the isotropic setting (2.2), (2.4) yields a second order convergence in s to the sharp interface 
limit; see e.g. ^52l |3l |45l [32l [28] for details. In particular, it needs to be assumed that the 
shape function ^, recall (3.7), satisfies ( 3.9a| b). Clearly, of our examples in (3.8) only the 
choices (iii) and (iv) satisfy this. 
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Table 8: Benehmark problem 2^^^ with ^ = 1 for PFEM. 
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80330 


22:12 mins 




512 


10" 


-4 


1.0521e-02 


1.3930e-02 


85218 


2:48 hours 




512 


10" 


-5 


4.1001e-03 


5.2299e-03 


85714 


27:46 hours 




512 


10" 


-6 


3.3949e-03 


6.4424e-03 


85842 


193:11 hours 



Table 9: Benchmark problem 



2(^) 



with ^ = 3 for PF°)^)^-FEM. 



28 





r 


r'' — r ll^oo 


\\u^ — u\\l<x> 


DOFs(r) 


CPU time 


32 


10" 


-2 


5.9808e-02 


3.4189e-02 


1213 


1 sees 


32 


10" 


-3 


4.2215e-02 


2.9898e-02 


1169 


3 sees 


32 


10" 


-4 


4.1433e-02 


2.9813e-02 


1217 


27 sees 


64 


10- 


-2 


2.7229e-02 


2.9493e-02 


2705 


1 sees 


64 


10- 


-3 


2.5199e-02 


1.7501e-02 


2505 


7 sees 


64 


10- 


-4 


2.6297e-02 


1.8341e-02 


2497 


1:06 mins 


128 


10- 


-2 


1.1307e-02 


3.0940e-02 


6133 


3 sees 


128 


10- 


-3 


1.0799e-02 


7.7496e-03 


5529 


21 sees 


128 


10- 


-4 


1.1358e-02 


8.4107e-03 


5489 


2:49 mins 


256 


10- 


-2 


5.9618e-03 


3.1694e-02 


15781 


8 sees 


256 


10- 


-3 


4.6057e-03 


3.1840e-03 


13165 


1:08 mins 


256 


10- 


-4 


4.9513e-03 


3.7938e-03 


12977 


10:08 mins 


512 


10- 


-2 


8.1252e-03 


3.2015e-02 


37069 


21 sees 


512 


10- 


-3 


1.9929e-03 


3.5308e-03 


35205 


3:44 mins 


512 


10- 


-4 


2.1894e-03 


1.7126e-03 


34149 


33:04 mins 


1024 


10- 


-2 


4.6688e-03 


3.2259e-02 


89301 


1:09 mins 


1024 


10- 


-3 


8.9176e-04 


3.7858e-03 


105993 


14:07 mins 


1024 


10- 


-4 


1.0202e-03 


7.8635e-04 


101517 


2:04 hours 



Table 10: Benchmark problem 2^^) for £ = 3 for PFEM 



2l/hf 


r 


\\r^ — r |ioo 


\\u^ — u\\l^ 


DOFs(r) 


CPU time 


32 


10- 


-2 


1.0373e-01 


4.0392e-02 


38661 


3:06 mins 


32 


10- 


-3 


1.2020e-01 


6.0675e-02 


38067 


23:07 mins 


32 


10- 


-4 


3.4579e-01 


1.0935e-01 


39826 


3:06 hours 


64 


10- 


-2 


2.4115e-02 


3.2858e-02 


162161 


25:31 mins 


64 


10- 


-3 


3.0456e-02 


1.9335e-02 


147787 


2:33 hours 


64 


10- 


-4 


4.7083e-02 


2.7765e-02 


148036 


17:10 hours 


128 


10- 


-2 


1.0819e-02 


3.1504e-02 


740635 


2:08 hours 


128 


10- 


-3 


1.1358e-02 


7.7776e-03 


603059 


12:43 hours 


128 


10- 


-4 


1.6654e-02 


1.0916e-02 


593187 


65:34 hours 


256 


10- 


-2 


8.2692e-03 


3.1338e-02 


3480836 


9:19 hours 


256 


10- 


-3 


5.6216e-03 


3.7231e-03 


2577960 


47:19 hours 


256 


10- 


-4 


6.4224e-03 


4.7100e-03 


2490788 


308:35 hours 



Table 11: Benehmark problem 3^^) with £ = 3 for PFEM. 



29 



£-1 




T 


1 1 XI 11-^ 


Ww'^ — 


DOFs(r) 


CPU time 


IGtt 


128 


10" 


-2 


3.6341e-02 


1.2977e-01 


15370 


13 sees 




128 


10" 


-3 


1.8831e-02 


1.8600e-02 


16226 


2:04 mins 




128 


10" 


-4 


1.3159e-02 


1.5797e-02 


15818 


14:12 mins 




128 


10" 


-5 


1.2840e-02 


1.5560e-02 


15890 


2:09 hours 




128 


10" 


-6 


1.2810e-02 


1.5510e-02 


15890 


19:13 hours 


32 TT 


256 


10" 


-2 


2.2563e-01 


2.1820e-01 


26362 


29 sees 




256 


10" 


-3 


6.1576e-03 


3.4619e-02 


34634 


5:04 mins 




256 


10" 


-4 


2.6943e-03 


7.5241e-03 


34938 


53:05 mins 




256 


10" 


-5 


2.8987e-03 


7.9607e-03 


34754 


5:23 hours 




256 


10" 


-6 


2.9172e-03 


1.1125e-02 


34818 


42:28 hours 


64 TT 


512 


10" 


-2 


3.8566e-01 


2.5703e-01 


58146 


1:27 mins 




512 


10" 


-3 


7.9018e-02 


1.0605e-01 


80010 


21:39 mins 




512 


10" 


-4 


1.0560e-02 


1.3817e-02 


85266 


2:51 hours 




512 


10" 


-5 


3.9822e-03 


5.3896e-03 


86002 


23:17 hours 




512 


10" 


-6 


3.1951e-03 


1.0395e-02 


86042 


154:02 hours 



Table 12: Benchmark problem 



4W 



with ^ = 3 for PFg-FEM. 



From now on we assume that (2.2) and (2.4) hold, and that /C+ = /C > 0. Then in 



place of (3.5a) and (3.6a) we consider 

'dwt + \ q{^) (/:^t = /C a It; , 

a Q{ip) w = £ {p + pis) (ft — ^ + ^\^) ^ 
in VLt-) where ^ : M ^ M is a second shape function that satisfies 

q{s)>Q V5G[-1,1], q{s) = q{-s) V^gM, J^g{y)dy = l. 



(4.6a) 
(4.6b) 



Here pi in (4.6b) is a correction term that is given by 

A a 



Pi = K 



(4.7) 



Pis) 



Q{y) ■ 



(4. 



where, see e.g. [45l Eq. (35)], 

K := Al - P($(5))] P($(5)) ds with 
Jr 

In (4.8) the function $ : M ^ M denotes the unique solution to 

^'\s) - ^'($(5)) = V 5 G M , lim $(5) = ±1 , [ s ^\s) ds = . 

For the above choices of ^, ^and pi, Almgren [3] formally showed second order convergence 
in the sense that the approximation of the Gibbs-Thomson law is of 0{e^)^ whereas in [45] 
it is formally established that 



30 



2l/hf 


r 


\\r^ — r||/,oo 




DOFs(r) 


CPU time 


32 


10" 


-2 


6.4835e-02 


3.5949e-02 


1225 


sees 


32 


10" 


-3 


4.3505e-02 


3.0646e-02 


1169 


5 sees 


32 


10" 


-4 


4.2440e-02 


3.0431e-02 


1217 


45 sees 


64 


10" 


-2 


2.9889e-02 


2.9534e-02 


2613 


1 sees 


64 


10" 


-3 


2.6450e-02 


1.8257e-02 


2521 


11 sees 


64 


10" 


-4 


2.7699e-02 


1.9103e-02 


2497 


1:40 mins 


128 


10" 


-2 


1.2103e-02 


3.1028e-02 


6109 


4 sees 


128 


10" 


-3 


1.1158e-02 


8.0117e-03 


5529 


27 sees 


128 


10" 


-4 


1.1813e-02 


8.7200e-03 


5489 


4:12 mins 


256 


10" 


-2 


6.4202e-03 


3.1800e-02 


15805 


10 sees 


256 


10" 


-3 


4.7724e-03 


3.3036e-03 


13221 


1:19 mins 


256 


10" 


-4 


5.1346e-03 


3.9258e-03 


13001 


12:16 mins 


512 


10" 


-2 


8.8182e-03 


3.2129e-02 


37161 


24 sees 


512 


10" 


-3 


2.0921e-03 


3.5487e-03 


35117 


4:19 mins 


512 


10" 


-4 


2.2635e-03 


1.7705e-03 


34165 


38:59 mins 


1024 


10" 


-2 


4.7558e-03 


3.2373e-02 


89269 


1:15 min 


1024 


10" 


-3 


9.5286e-04 


3.8105e-03 


105997 


16:11 mins 


1024 


10" 


-4 


1.0525e-03 


8.1355e-04 


101509 


2:25 hours 


2048 


10" 


-2 


4.6014e-03 


3.2830e-02 


295201 


6:13 min 


2048 


10" 


-3 


6.7747e-04 


3.8557e-03 


352153 


1:16 hours 


2048 


10" 


-4 


4.9986e-04 


3.6634e-04 


334785 


11:48 hours 



Table 13: Benchmark problem 4^^^ for ^ = 3 for PFEM. 



• the zero level set of the phase field function Lp approximates the interface V to 0{e^)^ 

• the temperature w in the phase field system approximates the temperature u in the 
sharp interface problem to 0{e^). 

In [32], for the special case ^{s) = |, the above second order approximation results are 
shown rigorously. In particular, on letting /C = a = 1, and on recalling that in their 
notation G{s) = c^P(5), it holds that the expression in [32l Eq. (1.6)] for the correction 
term pi is given by 

_ , ^ J^[G{1) - G{<^{s))]{l + $(.)) _ , ^ /Jl-P($(.))](l + $(.)) 
= \\ /[l-P($(s))](l + $(s))ds = Air, 



and so agrees exactly with (4.7). In addition, on assuming the stronger condition 



(4.9) 



31 



in place of ( 3.9a| b), the authors in [32] also show rigorously that the full phase field converges 
to second order. More precisely, in this case the first order correction to the phase field 
function is zero. Of course, the specific choice ^(5) = \ for the interpolation function in 
the equation for the temperature means that the overall phase field system considered in 
[32] is not thermodynamically consistent. We refer to [3l[32l|45] for the precise statements 
of these results. 



From now on we consider (4.6a,b) in the case that Q = which means that we return 
to (3.5a), (3.6a) in the isotropic case (2.2), (2.4) with /C+ = /C_. In particular, the phase 



field model ( 4.6a| b) is now thermodynamically consistent, i.e. it satisfies (3.12) with (2.2) 
(2.4) and /C+ = /C_. Since now ^= ^, it follows from (4.8) that 



K= / [l-P($(5))]P($(5))d5. 



In the case of the obstacle potential (3.4), so that $(5) = sin(5) for \s\ < | as in (3.18a) 
we get that 



K 



P(sin(s))] P(sin(s)) As. 



In particular, for the choice (3.8)(iii), when P(s) = (3 — 10 s"^ + 15 s + 8), we have that 



K = vr = 2"^^ 4817 n ^ 0.231 . 



65536 



(4.10) 



We refer to Table UM for computations for the benchmark problem 2^ with ^ = 1 for the 



phase field model (|4.6a|b) with the correction term (4.7) and (4. 



PFg-FEM the scheme PF^??f)-FEM but for the phase field model 



in (3.6a) replaced by (3.10) 



0\j . Here we d enote by 
d), (3.6a c) with p 



3.5a 



The numbers in Table [14] show that eight points across the interface are not enough to 
see convergence in the errors in practice. Here we recall that a similar conclusion can 

be drawn from the results reported in \2E, Table 1] , where a one-dimensional reformulation 
of a radially symmetric problem in is considered. With 16 points across the interface 
the error \\r^_ — r\\Loo in Table 14 appears to converge quadratically in ^, while only linear 



XI 



convergence can be seen in the error — t^||i,oo. We also note that the errors ||r^^ — rH^oo 



shown in Table [14] are significantly smaller than the corresponding errors in Tables [4] and 
[5] for the classical phase field model, i.e. for pi = 0, whereas the improvements in the 
temperature error — t^||i,oo are less pronounced. 



In the case of the quartic potential (3.3), so that $(5) = tanh(2 2 s) as in (3.18b), we 
get that 

K = V2 [[1- P(tanh(5))] P(tanh(5)) ds . 
Jr 



In particular, for the choice (3.8)(iii) we have that 



K=^^^ 0.352. 



840 



32 



£-1 




T 


1 1 XI 11-^ 


\\w^ — u\\l^ 


DOFs(T) 


CPU time 


47r 


32 


10" 


-4 


4.3481e-03 


2.7780e-02 


3394 


5:12 mins 




32 


10" 


-5 


1.9588e-03 


2.4734e-02 


3442 


39:41 mins 




32 


10" 


-6 


1.9743e-03 


2.4425e-02 


3474 


5:40 hours 




64 


10" 


-4 


1.5762e-03 


2.5180e-02 


11922 


24:35 mins 




64 


10" 


-5 


3.9787e-03 


2.1853e-02 


11906 


2:41 hours 




64 


10" 


-6 


4.5109e-03 


2.1533e-02 


11970 


29:05 hours 


V32 7r 


64 


10" 


-4 


9.3961e-03 


2.3676e-02 


8946 


19:06 mins 




64 


10" 


-5 


9.2783e-04 


1.7200e-02 


9074 


2:25 hours 




64 


10" 


-6 


1.1183e-03 


1.6723e-02 


9186 


19:37 hours 




128 


10" 


-4 


8.5330e-03 


2.2021e-02 


32226 


1:46 hours 




128 


10" 


-5 


9.0354e-04 


1.6249e-02 


32730 


13:02 hours 




128 


10" 


-6 


1.7709e-03 


1.5796e-02 


32586 


104:40 hours 


Stt 


64 


10" 


-4 


2.3635e-02 


2.9991e-02 


7074 


15:49 mins 




64 


10" 


-5 


6.8162e-03 


1.6316e-02 


7042 


1:46 hours 




64 


10" 


-6 


5.0485e-03 


1.5149e-02 


7066 


14:08 hours 




128 


10" 


-4 


1.9116e-02 


2.5065e-02 


24266 


1:12 hours 




128 


10" 


-5 


1.9668e-03 


1.2975e-02 


24746 


9:47 hours 




128 


10" 


-6 


4.8430e-04 


1.2136e-02 


24874 


65:32 hours 



Table 14: Benchmark problem 



2W 



with £ = 1 for PFjf)-FEM. 



while for the choice (3.8)(iv), on noting that P(s) 

K 



1(2 + 3 s - s^), it holds that 



_ 19\/2 
60 



0.448 . 



(4.11) 



We refer to Table 15 for computations for the benchmark problem 



phase field model (4.6a,b) with the correction term (4.7) and (4.1 



PFg-FEM the scheme PF^^^-FEM for the phase field model (|3.5a 



2^3 with ^ = 1 for the 
). H ere w e denote by 
-d), (3.6a-c) with p in 



(3.6a) replaced by (3.10), and with g given by (3.8)(iv). A computation for the scheme 



PF^^^-FEM with the implicit time discretization (3.20|) yielded very similar error numbers, 
and so we omit these results here. What can be 
again we have convergence of order at least 0{e^ 



clearly seen from Table [15] is that once 
in ||r^^ — tWloo. 



In order to visualize the performances of all the considered methods for the benchmark 
problem 2^^^ with ^ = 1, i.e. including the computations in Tables 14 and 15 for the second 



order accurate phase field model (4.6a,b), we present a log-log plot of the errors in the radius 



r against the necessary CPU time for all the entries in the appropriate tables in Figure [TO 
The plot appears to confirm that computations for the phase field model ( 4.6a[ b) are on 
average more efficient than computations for the standard phase field model, i.e. ( |4.6a[b) 
with pi = 0. However, the finer meshes needed for computations for (4.6a,b), recall Table |14t 
mean that due to CPU time constraints we cannot choose s as small as in the standard phase 
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£-1 




T 


1 1 XI 11-^ 


\\w^ — u\\l^ 


DOFs(T) 


CPU time 


47r 


32 


10" 


-4 


1.1945e-02 


4.1523e-02 


4066 


7:12 mins 




32 


10" 


-5 


1.7315e-02 


3.9842e-02 


4090 


52:56 mins 




32 


10" 


-6 


1.7860e-02 


3.9682e-02 


4082 


8:24 hours 




64 


10" 


-4 


1.1659e-02 


4.1790e-02 


14090 


29:48 mins 




64 


10" 


-5 


1.7071e-02 


3.9883e-02 


14130 


4:18 hours 




64 


10" 


-6 


1.7615e-02 


3.9686e-02 


14210 


43:32 hours 


V32 7r 


64 


10" 


-4 


5.1351e-03 


3.2940e-02 


10514 


18:46 mins 




64 


10" 


-5 


4.0145e-03 


2.9244e-02 


10578 


2:41 hours 




64 


10" 


-6 


5.0498e-03 


2.8980e-02 


10578 


25:57 hours 




128 


10" 


-4 


5.1726e-03 


3.2286e-02 


37194 


1:09 hours 




128 


10" 


-5 


3.7103e-03 


2.9080e-02 


37690 


10:37 hours 




128 


10" 


-6 


4.7916e-03 


2.8833e-02 


37626 


109:39 hours 


Stt 


64 


10" 


-4 


1.9976e-02 


3.2507e-02 


8082 


12:14 mins 




64 


10" 


-5 


3.4097e-03 


2.2891e-02 


8242 


2:12 hours 




64 


10" 


-6 


1.4844e-03 


2.2330e-02 


8306 


16:44 hours 




128 


10" 


-4 


1.9474e-02 


3.0950e-02 


27826 


1:15 hours 




128 


10" 


-5 


2.9987e-03 


2.2277e-02 


28298 


7:39 hours 




128 


10" 


-6 


1.0518e-03 


2.1768e-02 


28234 


65:58 hours 



Table 15: Benchmark problem 



2^ 



with 



1 for PFJ^^-FEM. 



field computations. Finally, the plot in Figure [10] once again underlines the superiority of 
the sharp interface algorithm PFEM over all the phase field methods. 



5 Numerical experiments for anisotropic problems 



In this section we present numerical simulations for the anisotropic Stefan problem ( 2.3a| -e). 
Here we always let = and ;5 = 7, where we recall (3.1). Moreover, we always 

choose A = a = 1 and, unless otherwise stated, we let 1C± = 1. In order to appreciate 



the computational challenges involved with the different experiments, we recall from (3.16) 



(3.17) and (3.11) that for accurate phase field calculations the following implications arise: 



~^ large 



e small 



/i/ , r small . 



Moreover, we note that for the fully anisotropic situation a formally second order accurate 



phase field model similarly to Section 4.3 involves a parameter piiV ^p) in place of (4.7) 



that depends on /3, 7 and on V(/:?, see ^52j,i53j. We stress that these approaches are not well 
analyzed so far, e.g. in the spirit of [45^, ^32]. In particular, to our knowledge there are no 
formal or rigorous results in the literature on the second order convergence of phase field 
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1e+06 



100000 




0.0001 



0,001 0.01 
Error in r 



Figure 10: (Benchmark problem 2^^^ 
LOO against the CPU time for t 

14 and 15 



Tables 
circles) . 



with = 1) Log-log scatter plot of ||r^^ — rH^oo and 
le entries in Tables|4|and|5|(PFj^^-FEM, blue rhombi), 



(PFg-FEM and PF^^^-FEM, green r 



(iv) 



lombi) and Table 



(PFEM, red 



models for the fully anisotropic Gibbs-Thomson law. Moreover, our numerical results in the 
isotropic case showed that the second order models do not give a large gain in computational 
efficiency. That is why all of our phase field computations in this section are for the standard 
phase field model ( 3.5a| ^d), ( 3.6a| -c). 



For the first simulations that we present we choose as anisotropy the regularized /^-norm 

2 

ANii: 7(p) = [5^ \p\^ + p] (1 - 5^)] ^ , with 5 = 0.3 . 

The radius of the initially circular seed Fq is chosen as 0.1, while we set i9 = 0, a = 0.03 and 
p = 0.01. The supercooling at the boundary = of Q = (—8, 8)^ is set to Ud = —2. 

Three numerical simulations for the scheme PF????-FEM with the interfacial parameter 



£~ — An can be seen in Figure |11[ It can be seen that varying the time discretization 
parameter r from 10~^ to 10~^ has a significant impact on the observed numerical results. 
However, the observed changes for the smallest value of r are small, which indicates that the 
simulation appears to be converging. Very similar results can be obtained for PF^^^-FEM 



with the implicit time discretization from (3.20), and so we omit them here 



We compare the above numerical experiments for the phase field method with three 



simulations for the sharp interface approximation PFEM in Figure 12, where we fix the 



spatial discretization parameters diS hr ^ ^hf = Here we observe that even for a 
very crude time discretization, the evolution is captured remarkably well, and there is very 
little variation in the numerical results from PFEM when r is decreased. Also note that it 
takes (less than) a second of CPU time with PFEM in order to get a good idea about the 
evolution of the growing crystal, while the phase field methods PF^^^-FEM and PF^^^-FEM 
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Figure 11: (PF^|?J-FEM, e'^ = 4 7r, ANii, i9 = 0, a = 0.03, p = 0.01, ud = -2, Q = 
(—8, 8)^) Snapshots of the solution at times t = 1, 2, 4, 5, 6. From top to bottom r = 10~^, 
k = 1 ^ 3. [These computations took 34 seconds, 7 minutes and 57 minutes, respectively.] 



take at least 400 times as long. 

For the next set of numerical experiments we use the hexagonal anisotropy 

3 

ANI2: -f{p) = ^l{R{^+ei)p), where l{p) = [pj + IQ-"" pj]'^ , (5.1) 
£=1 



and where R(9) = ( . . . ) denotes a clockwise rotation through the angle 9. Moreover, 

^ ^ sm cos J o o 7 

we use the parameters i9 = 1, a = 5 x 10~^, p = 0.01 and ud = —\ on the boundary 
= dO. of = (—2,2)^. The radius of the initially circular seed is again chosen as 
i?o = 0.1, and we set 







Ud 



1 - e^o 
Ud 



\x\ < Rq , 
— (1 - e^o-lo^l) i?o < |x| <H 

\x\> H , 



(5.2) 



with H := 2. 

Three numerical simulations for the scheme PF?j'?f-FEM with the interfacial parameter 



-1 _ 



16 TT can be seen in Figure 13, Observe that here we use a much smaller value of e 



than previously, because for larger values of e large mushy interfacial regions develop, which 
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Figure 12: (PFEM, ANii, ^ = 0, a = 0.03, p = 0.01, ud = -2, Q = (-8,8)^) Snapshots of 
the solution at times t = 1, 2, 4, 5, 6. From top to bottom r = 10~^, k = 1 ^ 3. [These 
computations took 1, 6 and 73 seconds, respectively.] 



means that the phase field simulations hold no value for the investigation of the underlying 
sharp interface problem. The creation of small localized mushy interfacial regions can be 



observed in Figure 13 for the run with r = 10~^, while the run with r = 10~^ shows larger 
such regions. In addition, in the latter run the phase field approximation of the growing 
crystal's surface reaches the external boundary dQ^ which results in the creation of artificial, 
nonphysical interfaces. Repeating these simulations for a smaller interfacial parameter s 
yields the results shown in Figure [T4| Now for sufficiently small values of the time step size 



r, the numerical results appear to be converging. 

We also repeat the last computation for the scheme PF^??J-FEM with the explicit time 
discretization from ( |3.22a ,b). Here any computation with a time step size r > 10~^ was 



unstable, and so in Figure |15| we only show a run for r = 10~ . We recall from Section |4.2 
see Tables |6] and [7[ that in the interest of accuracy uniform meshes should be employed for 
an explicit method. However, the large CPU times associated with a uniform grid mean 
that we are unable to complete the evolution within a reasonable amount of time. Hence 
in Figure 15 we use the same adaptive mesh strategy as in Figure 14 for the semi-implicit 
scheme PF^^?^^-FEM. Note that while the finest run in Figure agrees well with the results 
shown in Figure [T5| the very small time step size used for theTatter means that the explicit 
scheme takes about 40 times as long as the semi-implicit scheme to compute the evolution. 
Hence, without further code optimizations, the explicit scheme would need to be run in 
parallel on a cluster with at least 40 nodes to become competitive. 
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Figure 13: (PF^^J-FEM, e'^ = 16 ANI2, ^ = 1, a = 5 x 10-^ p = 0.01, ud = -|, 
Q = (—2,2)^) Snapshots of the solution at times t = 0.3, 0.4, 0.5, 0.6, 0.7. From top to 
bottom r = 10~^, = 2 ^ 4. [These computations took 29 seconds, 26 minutes and 381 
minutes, respectively.] 



for the scheme PFJ^^-FEM, 



When we repeat the simulations shown in Figures 13 andfl4 
then in the run for 

the boundary dQ. Hence we only present the simulations for s 
where we observe similar, but qualitatively quite different, results to the o nes shown in 
for the scheme PF^^^^^-FEM. In particular, the side arms in Figure 



16 TT large mushy interfacial regions appear, which quickly reach 

= 32 7r, see Figure 16 



14 



Figure 

be thinner than in Figure [141 and the convergence as r gets smaller appears 



16 



appear to 
be slower. 



We also repeat the computation from Figure 16 for the implicit time discretization from 
(3.20), see Figure 17, We observe that, in contrast to the conclusions that could be drawn 
from the isotropic experiments in Section [4| for the strongly anisotropic situation treated 
here there does not seem to be an advantage in using the implicit time discretization from 



(3.20) over the standard semi-implicit discretization ( 3.19a| b) from [16j 



In addition, we present three simulations for the same physical problem for the sharp 



interface approximation PFEM in Figure [T8| where we fix the spatial discretization pa- 
rameters as hr ^ hf = Here we observe once again that even for a very crude time 
discretization, the evolution is captured remarkably well, and there is very little variation in 
the numerical results from PFEM when r is decreased. We also draw particular attention 
to the differing CPU times between the sharp interface calculations in Figure [18] and the 



phase field simulations in Figures 14-17 
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Figure 14: (PF^^J-FEM, e'^ = 32^, ANI2, ^ = 1, a = 5 x 10-^ p = 0.01, ud = -|, 
Q = (—2,2)^) Snapshots of the solution at times t = 0.3, 0.4, 0.5, 0.6, 0.7. From top to 
bottom r = 10~^, A; = 2 ^ 4. [These computations took 38 seconds, 8 minutes and 69 
minutes, respectively.] 



In the remainder of this section we consider two simulations for the scheme PFEM, which 
on present computer hardware are virtually impossible to repeat to a desirable accuracy 
with the phase field method. The first experiment is with the physical parameters from [14[ 



Fig. 7], and so is for the one-sided quasi-stationary problem (2.3a-e) with i9 = /C_ = and 
with 7 as in (5.1). The remaining parameters are chosen as a = 10~^, p = 1.42 x 10~^ and 
Ud = —0.04 on the boundary = oiQ = (—4, 4)^. The radius of the initially circular 



seed is chosen as 0.05. See Figure [T9| for the results for different choices of the time step sizes 
r, and with hr ^ hf = fixed. We see that, as before, there is hardly any variation in the 
numerical results obtained from the three simulations with different values of r for PFEM. 
Moreover, we note that due the choice of the physical parameters much finer side branches 
appear in Figure [19] compared to the simulations in Figure [T8| To precisely capture these 



small structures within a phase field computation would require very small values for the 
interfacial parameter 6:, as well as correspondingly small discretization parameters hf and r; 
recall (3.17), and (3.11). Taken together this means that we are currently unable to present 



phase field computations for an evolution as shown in Figure [T9j 

The next computation is similar to the simulation shown in [12^ Fig. 14], where here we 
take as anisotropy 

ANI3: 7b) = {[9{P)T + [9{RiP)T + [9{R2p)Ty' , where g{p) := [pj + | {pj + pi 



(5.3) 
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Figure 15: (PF^^^^-FEM with the exphcit time discretization from (3.22a 



ANI2, i9 = 1, a = 5 X 10"^, p = 0.01, Ud = = (-2, 2f) Snapshots o 

times t = 0.3, 0.4, 0.5, 0.6, 0.7 for r = 10~^. [This computation took 44 hours. 



b), = 32 TT, 

" the solution at 





1 














0^ and i?2 := 


1 


0). The 


^ 







-1 


0^ 


the 


ri^ 


iht of Figure 3 in 


[12J. 


Moreover, 



we use the parameters i9 = l,a = 10 ^,p = 0.01 and Ud = on the boundary 
ofn= (-4,4)3 



10~^, p = 0.01 and = — ^ on the boundary = dfl 
The radius of the initiaUy spherical seed is chosen as i?o = 0.1, and we let 
Three simulations for these parameters, with the spatial 
are presented in Figure We observe 



Uq be defined by (5.2) with H = 4. 
discretization parameters fixed as /ir ^ 5/i/ 

that the three simulations all show the same general shape of the growing six-armed crystal, 
and for the smallest value of r the results appear to have converged. We also note that the 
small oscillations in the solution for the simulation with the largest time step size disappear 
as T is decreased. In order to demonstrate the good mesh properties of the parametric 

two details of the triangulated approximation of r{t) 

We recall 
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method PFEM, we show in Figure 

at times t = 0.2 and t = 0.3 for the finest time discretization in Figure [20 



that the algorithm PFEM does not employ any mesh-redistribution or mesh-smoothing 
methods. Rather it relies solely on local mesh refinements, where individual elements of the 
triangulation become too large, see [12' §5.2] for more details. The quality of the meshes 



shown in Figure 21 is excellent. 



We recall that a simulation for the phase field algorithm PF^^J-FEM for the same physical 
parameters has recently been performed in [16^ Fig. 23]. As a comparison to the sharp 
interface calculations from Figure 
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we present the results for the scheme PF^^J-FEM in 
Figure [22j We note that the evolution shown for the phase field approximation in Figure [22 
is qualitatively very different from the sharp interface simulations in Figure [20) In all 
likelihood the physically challenging parameters for the computation in Figure 22 mean 
that, both in terms of the discretization parameters for the given = 16 tt, e.g. the time 
step size r, as well as in terms of the interfacial parameter e itself, the shown numerical 
results are still far from the true underlying solution to the sharp interface problem ( 2.3a| -e). 
Of course, a detailed numerical study into this question is not yet possible due to the long 
time that such computations would take. 
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Figure 16: (PFj^-FEM, s'^ = 32 tt, ANI2, i9 = 1, a = 5 x 10-^ p = 0.01, = 

= (—2,2)^) Snapshots of the solution at times t = 0.3, 0.4, 0.5, 0.6, 0.7. From top to 
bottom r = 10~^, /c = 2 ^ 4. [These computations took 34 seconds, 16 minutes and 101 
minutes, respectively.] 



Conclusions 

While numerical simulations for phase field models in general show qualitatively correct 
behaviour, often such numerical results are far away from the true sharp interface evolution. 
In order to obtain accurate simulations, the interface width 6:, as well as the spatial and 
temporal discretization parameters need to be chosen sufficiently small. However, reducing 
these parameters to reach an acceptable accuracy often requires very large computing times 
on even the most advanced of today's desktop computers. 

Direct sharp interface approximations, on the other hand, can provide a computationally 
cheap method to compute interface evolutions in materials science accurately. An example 
of such an algorithm is PFEM from [T2l [T3] . In the computations presented in this paper we 
have seen that even for very crude discretization parameters, the algorithm PFEM provides 
surprisingly accurate approximations. The computational time needed to compute these 
sharp interface approximations is often negligible compared to the CPU times necessary for 
a corresponding phase field simulation. 

The main problem of phase field methods is that the asymptotic error in 6:, which in 
general is not known, plays a significant role in determining the accuracy of phase field sim- 
ulations. Small values of 6:, in turn, require very small discretization parameters. Similarly, 
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Figure 17: (PF^^^^-FEM with the imphcit time discretization from (|3.20|), e 



-- 32 TT, ANI2, 

1, a = 5 X 10~^, p = 0.01, Ud = — |, ^ = (—2, 2)^) Snapshots of the solution at times 
t = 0.3, 0.4, 0.5, 0.6, 0.7. From top to bottom r = 10"^, = 2 ^ 4. [These computations 
took 43 seconds, 23 minutes and 119 minutes, respectively.] 



in second order convergent isotropic phase field models with a correction term, where the 
asymptotic error in s may be assumed to be relatively smaller than in classical phase field 
models, very small discretization parameters need to be employed in order to benefit from 
the smaller asymptotic error in practice. All of these issues do not arise in sharp interface 
approximations. 

The main advantage of phase field methods over direct front tracking methods is that 
they intrinsically allow for topological changes. However, for the problem of solidification 
and dendritic growth as considered in this paper, topological changes in general do not 
occur during the simulation of physically relevant evolutions for single crystals. 

In the past, researchers and scientists may have been discouraged from applying front 
tracking methods because of the diflSculties in implementing such methods and because of 
the deterioration of the mesh quality as the approximated sharp interface evolves in time. 
However, assembling the system matrices in parametric finite element methods for evolving 
manifolds is not much different from the assembly in standard Cartesian problems, see e.g. 
[38^ 39J. Of course, the coupling between a lower dimensional parametric mesh and a bulk 
mesh is nontrivial, but successful implementations have been used in e.g. [Ml El [13] . 
Moreover, the good mesh properties of the scheme PFEM from [121 113 i^ean that a good 
mesh quality is maintained throughout the numerical simulations, and no remeshing is 
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Figure 18: (PFEM, ANI2, i9 = 1, a = 5 x 10-^ p = 0.01, ud = = (-2, 2)^) Snapshots 

of the solution at times t = 0.3, 0.4, 0.5, 0.6, 0.7. From top to bottom r = 10"^, /c = 2 ^ 4. 
[These computations took 1, 9 seconds and 87 seconds, respectively.] 

required in practice. In fact, all the simulations presented in this paper were performed 
without any remeshing, see [12j for more details. 

We can summarize our main conclusions as follows: 

(CI) The parametric front tracking method PFEM is more accurate and computationally 
more efficient than phase field methods. 

(C2) For isotropic problems, implicit time discretizations for phase field models are often 
more accurate than semi-implicit time discretizations. 

(C3) Explicit time discretizations for phase field models need very small time steps in 
practice, and hence computations with explicit schemes are only competitive if run in 
parallel on large clusters. 

(C4) Second order accurate phase field models need finer discretization parameters than 
classical phase field models in order to demonstrate their superior approximation 
properties in practice. 

Finally we note that while the focus of this paper has been the problem of dendritic 
solidification, it is to be expected that similar conclusions can be drawn when considering the 
respective merits of phase field models and sharp interface methods for other free boundary 
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Figure 19: (PFEM, ANI2, ^ = 0, a = 10-^ p = 1.42 x 10-^ ud = -0.04, Q = (-4,4)^) 
Snapshots of the solution at times t = 5, 10, 20, 30, 40. From top to bottom r = 10~^, 
k = 1 ^ 3. [These computations took 6, 124 and 781 minutes, respectively.] 

problems in materials science, physics and biology. As possible examples we refer to epitaxial 
growth, surface diffusion, thermal grooving, sintering, vesicle dynamics and two phase flow. 

It is our hope that the comparisons presented in this paper encourage a discussion about 
the merits of phase fleld methods in general, and of the possible advantages of using sharp 
interface approximations instead. 
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